#include<stdio.h>
#include<stdlib.h>
#include<string.h> 
#include<math.h>
#include<time.h>

//宏定义
#define MAX_POINTS 1000  
#define MAX_ITERATIONS 10000  
#define LEARNING_RATE 0.0001  
#define EPSILON 1e-6  

// 不同的函数类型
typedef enum {
    LINEAR = 1,      
    QUADRATIC = 2,   
    CUBIC = 3        
} FunctionType;

// 两种的算法
typedef enum {
    LS,       // 最小二乘法
    GD        // 梯度下降法
} AlgorithmType;

// 步骤一：生成数据集
int generate_data(const char* filename, FunctionType type, double a, double b, double c, double d) {
    FILE* fp = fopen(filename, "w");
    if (!fp) { 
        perror("打开文件失败"); 
        return 0; 
    }
    
    srand(time(NULL));
    fprintf(fp, "x,y\n"); 
    
    for (int i=0; i<MAX_POINTS; i++) {
        //判断x的取值范围
        double x = (double)(rand()%20000-10000)/100; 
        double y;
        switch(type) {
            case LINEAR:   y = a*x + b; break;
            case QUADRATIC: y = a*x*x + b*x + c; break;
            case CUBIC:    y = a*x*x*x + b*x*x + c*x + d; break;
            default: 
                fprintf(stderr, "无效函数类型\n"); 
                fclose(fp); 
                return 0;
        }
        fprintf(fp, "%.6f,%.6f\n", x, y); 
    }
    fclose(fp);
    printf("数据集已生成至 %s\n", filename);
    return 1;
}

// 步骤二：数据集的加载
int load_data(const char* filename, double* x, double* y) {
    FILE* fp = fopen(filename, "r");
    if (!fp) { 
        perror("加载文件失败"); 
        return 0; 
    }
    
    char line[100];
    fgets(line, sizeof(line), fp); 
    int count=0;
    
    while (fgets(line, sizeof(line), fp) && count<MAX_POINTS) {
        if (sscanf(line, "%lf,%lf", &x[count], &y[count]) == 2) count++;
    }
    fclose(fp);
    return count;
}

// 步骤三（1）：最小二乘法拟合
void least_squares(int n, int m, double* x, double* y, double* params) {
    double** A = (double**)malloc((m+1)*sizeof(double*));
    double* B = (double*)malloc((m+1)*sizeof(double));
    
    for (int i=0; i<=m; i++) {
        A[i] = (double*)malloc((m+1)*sizeof(double));
        for (int j=0; j<=m; j++) A[i][j] = 0;
        B[i] = 0;
    }
    
    for (int i=0; i<n; i++) {
        double xi = x[i], yi = y[i];
        for (int j=0; j<=m; j++) {
            double xi_j = pow(xi, j);
            B[j] += yi * xi_j;
            for (int k=0; k<=m; k++) {
                A[j][k] += xi_j * pow(xi, k);
            }
        }
    }
    
    // 高斯消元法
    for (int i=0; i<=m; i++) {
        int max_row = i;
        for (int j=i+1; j<=m; j++) {
            if (fabs(A[j][i]) > fabs(A[max_row][i])) max_row = j;
        }
     
        double* tmp = A[i]; A[i] = A[max_row]; A[max_row] = tmp;
        double t = B[i]; B[i] = B[max_row]; B[max_row] = t;
        
   
        double pivot = A[i][i];
        for (int j=i; j<=m; j++) A[i][j] /= pivot;
        B[i] /= pivot;
        
     
        for (int j=0; j<=m; j++) {
            if (j != i && fabs(A[j][i]) > EPSILON) {
                double factor = A[j][i];
                for (int k=i; k<=m; k++) A[j][k] -= factor * A[i][k];
                B[j] -= factor * B[i];
            }
        }
    }
    
  
    for (int i=0; i<=m; i++) params[i] = B[i];
    
 
    for (int i=0; i<=m; i++) free(A[i]);
    free(A); free(B);
}

// 步骤三（2）：梯度下降法拟合
void gradient_descent(int n, int m, double* x, double* y, double* params) {
    for (int i=0; i<=m; i++) params[i] = 0; // 初始化参数
    
    for (int iter=0; iter<MAX_ITERATIONS; iter++) {
        double *grad = (double *)malloc((m + 1) * sizeof(double));
        
        for (int i = 0; i < m + 1; i++) {
            grad[i] = 0;
        }
    
        for (int i=0; i<n; i++) {
            double xi = x[i], yi = y[i], pred=0;
            for (int j=0; j<=m; j++) pred += params[j] * pow(xi, j);
            double error = pred - yi;
            for (int j=0; j<=m; j++) grad[j] += error * pow(xi, j);
        }
        
        
        for (int j=0; j<=m; j++) {
            params[j] -= LEARNING_RATE * (2.0/n) * grad[j];
        }
        
     
        free(grad);
        
      
        double norm = 0;
        for (int j=0; j<=m; j++) norm += grad[j]*grad[j];
        if (sqrt(norm) < EPSILON) {
            printf("梯度下降在%d次迭代后收敛\n", iter+1);
            break;
        }
    }
}

// 步骤四：计算拟合误差
double calculate_sse(int n, int m, double* x, double* y, double* params) {
    double sse=0;
    for (int i=0; i<n; i++) {
        double pred=0;
        for (int j=0; j<=m; j++) pred += params[j] * pow(x[i], j);
        sse += pow(pred - y[i], 2);
    }
    return sse;
}

// 步骤五：计算决定系数
double calculate_r_squared(int n, int m, double* x, double* y, double* params) {
    double sse = calculate_sse(n, m, x, y, params);
    double mean_y = 0;
    for (int i=0; i<n; i++) mean_y += y[i];
    mean_y /= n;
    
    double sst = 0;
    for (int i=0; i<n; i++) {
        sst += pow(y[i] - mean_y, 2);
    }
    
    return 1.0 - (sse / sst);
}

// 步骤六：结果输出
void print_results(AlgorithmType alg, FunctionType type, double* params, double sse, double r_squared) {
    printf("\n=== 拟合结果 ===\n");
    printf("算法: %s\n", alg==LS?"最小二乘法":"梯度下降法");
    
    printf("函数形式: ");
    switch(type) {
        case LINEAR:
            printf("y = %.6fx + %.6f\n", params[1], params[0]);
            break;
        case QUADRATIC:
            printf("y = %.6fx2 + %.6fx + %.6f\n", params[2], params[1], params[0]);
            break;
        case CUBIC:
            printf("y = %.6fx3 + %.6fx2 + %.6fx + %.6f\n", params[3], params[2], params[1], params[0]);
            break;
    }
    
    // 输出评估指标
    printf("SSE (残差平方和): %.6f\n", sse);
    printf("R2 (决定系数): %.6f\n", r_squared);
    
    // 评估拟合质量
    if (r_squared > 0.95) {
        printf("拟合质量: 优秀 (R2 > 0.95)\n");
    } else if (r_squared > 0.8) {
        printf("拟合质量: 良好 (0.8 < R2 ≤ 0.95)\n");
    } else if (r_squared > 0.6) {
        printf("拟合质量: 一般 (0.6 < R2 ≤ 0.8)\n");
    } else {
        printf("拟合质量: 较差 (R2 ≤ 0.6)\n");
    }
}

int main(int argc, char* argv[]) {
    if (argc < 2) { // 至少需要一个命令参数
        printf("用法:\n");
        printf("生成数据: %s gendat <文件> <类型(1-3)> <a> <b> <c> <d>\n", argv[0]);
        printf("拟合数据: %s fitting <文件> <算法(ls/gd)> <类型(1-3)>\n", argv[0]);
        return 1;
    }
    
    if (strcmp(argv[1], "gendat") == 0) { 
        if (argc < 8) {
            printf("错误: 生成数据需要至少7个参数\n");
            printf("用法: %s gendat <文件> <类型(1-3)> <a> <b> <c> <d>\n", argv[0]);
            return 1;
        }
        
        int type_val = atoi(argv[3]);
        if (type_val < 1 || type_val > 3) {
            printf("错误: 函数类型必须是1(线性)、2(二次)或3(三次)\n");
            return 1;
        }
        FunctionType type = (FunctionType)type_val;
        
        double a = atof(argv[4]);
        double b = atof(argv[5]);
        double c = atof(argv[6]);
        double d = atof(argv[7]);
        
        if (!generate_data(argv[2], type, a, b, c, d)) {
            printf("数据生成失败\n");
            return 1;
        }
    } 
    else if (strcmp(argv[1], "fitting") == 0) { 
        if (argc < 5) {
            printf("错误: 拟合数据需要至少5个参数\n");
            printf("用法: %s fitting <文件> <算法(ls/gd)> <类型(1-3)>\n", argv[0]);
            return 1;
        }
        
        double x[MAX_POINTS], y[MAX_POINTS];
        int n = load_data(argv[2], x, y);
        if (n==0) {
            printf("无法加载数据，程序退出\n");
            return 1;
        }
        
        AlgorithmType alg;
             if      (strcmp(argv[3], "ls")==0) alg=LS;
        else if (strcmp(argv[3], "gd")==0) alg=GD;
        else { 
            printf("错误: 无效算法，必须是ls(最小二乘法)或gd(梯度下降法)\n"); 
            return 1; 
        }
        
        int type_val = atoi(argv[4]);
        if (type_val < 1 || type_val > 3) {
            printf("错误: 函数类型必须是1(线性)、2(二次)或3(三次)\n");
            return 1;
        }
        FunctionType type = (FunctionType)type_val;
        
        int m = type-1; 
        double params[4] = {0}; 
        
        printf("开始拟合...\n");
        switch(alg) {
            case LS:
                least_squares(n, m, x, y, params);
                break;
            case GD:
                gradient_descent(n, m, x, y, params);
                break;
        }
        
        double sse = calculate_sse(n, m, x, y, params);
        double r_squared = calculate_r_squared(n, m, x, y, params);
        print_results(alg, type, params, sse, r_squared);
    } 
    else {
        printf("错误: 未知命令\n");
        return 1;
    }
    
    return 0;
}
